# AN EXTRAVERSION AFFAIR: AFRICAN STATES AND LARGESCALE LAND ACQUISITION SUBNATIONAL LOCATIONS
# Version of the script: 12/04/2022
# C�cile Richetta, University of Geneva, Department of Political Science and International Relations 

# Figures in main text
  # Figure 1(a) 
  # Figure 1(b) was done in ArcGIS-Pro, not replicable. 
  # Figure 2
  # Figure 3
  # Figure 4(a) to Figure 4(e)


options(scipen=999)
options(digits=5)

library(base)
library(data.table)
library(foreign)
library(plyr)
library(dplyr)
library(Matrix)
library(splm)
library(spdep)
library(ggplot2)
library(arm)
library(dotwhisker)
library(hrbrthemes)
library(lubridate)

setwd("")
# data on land deals 
load("data.RData")
# weight matrix
W <- readMM(file = "W.mtx") # matrix
W <- as(W, "dgCMatrix") # to dgc format
W_lw <- mat2listw(W, style="W")

# Figure 1(a): number of lsla, new and cumulated, by year
  data2 <- data[, c(2, 69)]
  dataf <- aggregate(data2$newlla, by=list(Year=data2$year), FUN=sum)
  names(dataf)[2] <- "avgNEW"
  dataf$v1 <- as.character(as.numeric(dataf$Year))
  dataf$v2 <- paste(dataf$v1, "01-01", sep="-")
  dataf$v2 <- ifelse(dataf$Year==2018, "2018-12-31", dataf$v2)
  dataf$Year <- ymd(dataf$v2)
  
  # A few constants
  newColor <- "grey12"
  sec_breaks <- c(0, 120, 240, 360, 480)
  
  ggplot(dataf, aes(x=Year)) +
    geom_line( aes(y=avgNEW), size=1, color=newColor) +
    scale_y_continuous(
      # Features of the first axis
      name = "Newly Signed LSLA",
    ) + 
    theme_ipsum() +
    theme(
      axis.title.y = element_text(color = newColor, size=11),
      axis.title.y.right = element_text(color = newColor, size=11),
      plot.margin = unit(x = c(-1, 1, 0, 0), units = "mm"),
      axis.text.x = element_text(size=9)
    ) +
    ggtitle("")
  ggsave("figure1a.jpeg", width = 8, height = 10, units = "cm", dpi=1200)

# Figure 1(b): map of LSLA in units of analysis done in ArcGIS-Pro, script not available. 
  data2 <- data[, c(1, 69)]
  datald <- aggregate(data2$newlla, by=list(Year=data2$ID), FUN=sum)
  datald$x <- ifelse(datald$x>1, 1, datald$x)
  write.csv(datald, "unit-affected.csv", row.names = F)
  # the rest of the figure was done in ArcGIS : the dataset estimated in the lines above was merged with the shapefile unit IDs
  
# Figure 2: results of SPMLE analysis 
  setwd("") # directory where model objects are saved
  Model1A <- readRDS("Model1A.rds")
  Model1B <- readRDS("Model1B.rds")
  
  dwplot(list(Model1B, Model1A),
         vline = geom_vline(
           xintercept = 0,
           colour = "grey60",
           linetype = 2
         ),
         # plot line at zero _behind_ coefs
         dot_args = aes(size=1.3),
         whisker_args = aes(size=0.8),
         vars_order = c("nle_l1", "pop_l1", "urban_l1", "crop_l1", "rainsum_l1", "water_l1", "polvio_l1", "rho", "gamma")
  ) %>% # plot line at zero _behind_coefs
    relabel_predictors(
      c(
        nle_l1 = "Infrastructural Capacity",
        pop_l1 = "Population",
        urban_l1 = "Urbanization",
        crop_l1 = "Cropland",
        rainsum_l1 = "Rainfall",
        water_l1 = "Water",
        polvio_l1 = "Instability",
        rho = "Rho",
        gamma = "Gamma"
      )
    ) + # plot line at zero _behind_ coefs
    theme_bw(base_size = 20) + xlab("Coefficient Estimates") + ylab("") +
    theme(
      plot.margin = unit(x = c(1, 2, 1, -7), units = "mm"),
      legend.position="top",
      legend.box.margin = unit(x = c(0,0,-5,0), units = "mm"),
      legend.title = element_blank(),
      legend.margin=margin(t = 0, unit='cm'),
      axis.text.x = element_text(size=10),
      axis.title.x = element_text(size=11),
      axis.text.y = element_text(size=10),
      legend.text = element_text(size = 8)
    )  + 
    scale_colour_grey(
      start = .3,
      end = .7,
      labels = c("Model 1B", "Model 1A")
    ) 
  ggsave("figure2.jpeg", width = 16, height = 7, units = "cm", dpi=1200)

#Figure 3(a)
  Model2 <- readRDS("Model2.rds") # not significant but coherent
  impac1 <- impacts(Model2, listw = mat2listw(W, style = "W"), time = 11)
  a <- summary(impac1, zstats=TRUE, short=TRUE)
  
  x <- c("model", "term", "estimate", "std.error")
  varnames<-c("IC X Population", "IC X Urbanization", "IC X Cropland", "IC X Rainfall", "IC X Water", "IC X Instability")
  varnames2<-c("Infrastructural Capacity (IC)", "Population", "Urbanization", "Cropland",  "Rainfall", "Water", "Instability")
  
  direct1 <- data.frame(matrix(ncol = 4, nrow=7))
  colnames(direct1) <- x
  direct1$estimate <- summary(impac1)$direct$statistics[(1:7),1]
  direct1$std.error <- summary(impac1)$direct$statistics[(1:7),2]
  direct1$term <- varnames2
  direct1$model <- 0
  
  direct2 <- data.frame(matrix(ncol = 4, nrow=6))
  colnames(direct2) <- x
  direct2$estimate <- summary(impac1)$direct$statistics[(50:55),1]
  direct2$std.error <- summary(impac1)$direct$statistics[(50:55),2]
  direct2$term <- varnames
  direct2$model <- 0
  
  indirect1 <- data.frame(matrix(ncol = 4, nrow=7))
  colnames(indirect1) <- x
  indirect1$estimate <- summary(impac1)$indirect$statistics[(1:7),1]
  indirect1$std.error <- summary(impac1)$indirect$statistics[(1:7),2]
  indirect1$term <- varnames2
  indirect1$model <- 1
  
  indirect2 <- data.frame(matrix(ncol = 4, nrow=6))
  colnames(indirect2) <- x
  indirect2$estimate <- summary(impac1)$indirect$statistics[(50:55),1]
  indirect2$std.error <- summary(impac1)$indirect$statistics[(50:55),2]
  indirect2$term <- varnames
  indirect2$model <- 1
  
  total1 <- data.frame(matrix(ncol = 4, nrow=7))
  colnames(total1) <- x
  total1$estimate <- summary(impac1)$total$statistics[(1:7),1]
  total1$std.error <- summary(impac1)$total$statistics[(1:7),2]
  total1$term <- varnames2
  total1$model <- 2
  
  total2 <- data.frame(matrix(ncol = 4, nrow=6))
  colnames(total2) <- x
  total2$estimate <- summary(impac1)$total$statistics[(50:55),1]
  total2$std.error <- summary(impac1)$total$statistics[(50:55),2]
  total2$term <- varnames
  total2$model <- 2
  
  df <- rbind( direct1, direct2,  indirect1, indirect2, total1, total2)
  df <- data.table(df)

  dwplot(df,
         vline = geom_vline(
           xintercept = 0,
           colour = "grey60",
           linetype = 2
         ),
         # plot line at zero _behind_ coefs
         dot_args = aes(size=1.3),
         whisker_args = aes(size=0.8),
         vars_order = c( "Infrastructural Capacity (IC)", "Population", "Urbanization", 
                        "Cropland", "Rainfall", "Water", "Instability", 
                        "IC X Population", "IC X Urbanization", "IC X Cropland",  "IC X Rainfall", "IC X Water",
                         "IC X Instability")
  ) + # plot line at zero _behind_ coefs
    theme_bw(base_size = 20) + xlab("Impact Estimates") + ylab("") +
    theme(
      plot.margin = unit(x = c(1, 1, 1, -5), units = "mm"),
      legend.position="top",
      legend.box.margin = unit(x = c(0,0,-5,0), units = "mm"),
      legend.title = element_blank(),
      legend.margin=margin(t = 0, unit='cm'),
      axis.text.x = element_text(size=10),
      axis.title.x = element_text(size=11),
      axis.text.y = element_text(size=10),
      legend.text = element_text(size = 8)
    )  + 
    scale_colour_grey(
      start = .3,
      end = .7,
      labels = c("Total", "Indirect", "Direct")
    ) 
  ggsave("Figure3.jpeg", width = 16, height = 10, units = "cm", dpi=1200)
  
# Figure 4(a) to (f): marginal effects with meplot()
  meplot <- function(model,var1,var2,int,vcov,ci=.95,
                     xlab=var2,ylab=paste("Marginal Effect of",var1),
                     main="Marginal Effect Plot",
                     me_lty=1,me_lwd=1,me_col="black",
                     ci_lty=1,ci_lwd=.5,ci_col="black",
                     yint_lty=2,yint_lwd=1,yint_col="black"){
    require(ggplot2)
    alpha <- 1-ci
    z <- qnorm(1-alpha/2)
    beta.hat <- coef(model)
    cov <- vcov
    z0 <- seq(min(model.frame(model)[,var2],na.rm=T),max(model.frame(model)[,var2],na.rm=T),length.out=1000)
    dy.dx <- beta.hat[var1] + beta.hat[int]*z0
    se.dy.dx <- sqrt(cov[var1,var1] + z0^2*cov[nrow(cov),ncol(cov)] + 2*z0*cov[var1,ncol(cov)])
    upr <- dy.dx + z*se.dy.dx
    lwr <- dy.dx - z*se.dy.dx
    ggplot(data=NULL,aes(x=z0, y=dy.dx)) +
      labs(x=xlab,y=ylab,title=main) +
      geom_line(aes(z0, dy.dx),size = me_lwd, 
                linetype = me_lty, 
                color = me_col) +
      geom_line(aes(z0, lwr), size = ci_lwd, 
                linetype = ci_lty, 
                color = ci_col) +
      geom_line(aes(z0, upr), size = ci_lwd, 
                linetype = ci_lty, 
                color = ci_col) +
      geom_hline(yintercept=0,linetype=yint_lty,
                 size=yint_lwd,
                 color=yint_col)
  }
  
  meplot(model=Model2,var1="nle_l1",var2="pop_l1",int="nle_l1:pop_l1",vcov=vcov(Model2),
         xlab="Population",ylab="Marginal Effect IC",
         main="") + 
    theme_classic()                        # Use the classic theme.
  ggsave("Figure4a.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)
  
  meplot(model=Model2,var1="nle_l1",var2="urban_l1",int="nle_l1:urban_l1",vcov=vcov(Model2),
         xlab="Urbanization",ylab="",
         main="") + 
    theme_classic() +
    theme(
      plot.margin = unit(x = c(1, 3, 1, -3), units = "mm"))
  ggsave("Figure4b.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)
  
  meplot(model=Model2,var1="nle_l1",var2="crop_l1",int="nle_l1:crop_l1",vcov=vcov(Model2),
         xlab="Cropland",ylab="",
         main="") + 
    theme_classic() 
  ggsave("Figure4c.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)  
  
meplot(model=Model2,var1="nle_l1",var2="rainsum_l1",int="nle_l1:rainsum_l1",vcov=vcov(Model2),
         xlab="Rainfall",ylab="Marginal Effect IC",
         main="") + 
    theme_classic()
  ggsave("Figure4d.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)
  
  meplot(model=Model2,var1="nle_l1",var2="water_l1",int="nle_l1:water_l1",vcov=vcov(Model2),
         xlab="Water",ylab="",
         main="") + 
    theme_classic()
  ggsave("Figure4e.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)
  
  meplot(model=Model2,var1="nle_l1",var2="polvio_l1",int="nle_l1:polvio_l1",vcov=vcov(Model2),
         xlab="Instability",ylab="",
         main="") + 
    theme_classic() 
  ggsave("Figure4f.jpeg", 
         width = 5, height = 5, units = "cm", dpi=1200)

  
  
  
  # ---------------------------------------- END SCRIPT ------------------------------------------#